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DISCRETE  FILTERING  ON  UNSTRUCTURED  GRIDS  BASED 
ON  LEAST-SQUARES  GRADIENT  RECONSTRUCTION 


A.  HASELBACHER 

Center  for  Simulation  of  Advanced  Rockets 
University  of  Illinois  at  Urbana- Champaign 
Urbana,  IL  61801 


1.  Introduction 

The  equations  for  Large-Eddy  Simulation  (LES)  of  turbulent  flows  are  for- 
mally derived  by  applying  a low-pass  filter  to  the  Navier-Stokes  equations. 
In  doing  so,  it  is  often  tacitly  assumed  that  the  filtering  and  differentiation 
operations  commute.  This  assumption  is  invalid  if  the  filter  width  is  not 
uniform — as  is  the  case  if  wall-bounded  flows  are  computed — unless  special 
filter  operators  are  constructed,  see,  e.g,  Vasilyev  et  al.  (1998). 

Recent  work  by  Marsden  et  al.  (2000)  resulted  in  a framework  for  the 
construction  of  filters  on  unstructured  grids  which  commute  with  differ- 
entiation to  a potentially  arbitrarily  high  order.  They  also  demonstrated 
a filter  operator  with  a second-order  commutation  error.  However,  their 
method  appears  to  be  quite  complicated  in  its  construction,  particularly 
in  three  dimensions  and  near  boundaries.  Furthermore,  it  is  dependent  on 
geometric  comparisons  with  user-specified  parameters. 

The  goal  of  the  present  work  is  to  develop  a simpler  filtering  method 
than  that  of  Marsden  et  al.  (2000).  The  new  filtering  method  is  based  on  the 
following  observation:  The  conditions  for  filtering  a function  to  a given  order 
of  commutation  error  derived  by  Vasilyev  et  al.  (1998)  are  formally  identi- 
cal to  the  conditions  for  reconstructing  the  gradient  of  a function  to  a given 
order  of  truncation  error.  In  other  words,  the  construction  of  filtering  oper- 
ators may  be  reinterpreted  as  the  construction  of — suitably  reformulated — 
gradient-reconstruction  methods.  This  apparently  trivial  observation  has 
important  consequences  because  the  reconstruction  of  gradients  is  central 
to  many  flow-solution  methods  on  unstructured  grids  and  is  well  under- 
stood, see  Haselbacher  and  Blazek  (2000). 
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2.  Least-Squares  Gradient  Reconstruction 

The  least-squares  gradient-reconstruction  procedure  originally  developed 
by  Barth  (1991)  is  based  on  approximating  the  variation  along  an  edge 
linking  vertices  0 and  i by  a truncated  Taylor  series,  e.g.,  for  a linear  ap- 
proximation, 

4>i  = <t>o  + (V<£)o  • Ar0i,  (1) 

where  Aroj  = r*  — ro  and  r = {x,y}1.  The  application  of  Eq.  (1),  or 
corresponding  higher-order  approximations,  to  all  edges  incident  to  vertex 
0 gives  a system  of  linear  equations  for  the  derivatives  at  vertex  0, 

Ax  = b,  (2) 


where  A is  a do  x no  matrix  of  geometrical  terms,  x is  an  no-vector  contain- 
ing derivatives,  and  b is  a do-vector  of  function  values,  with  no  being  the 
number  of  derivatives  reconstructed  and  do  denoting  the  degree  of  vertex 
0.  Since  there  are  usually  more  incident  edges  than  derivatives,  Eq.  (2)  is 
solved  for  x in  a least-squares  fashion. 

A general  closed-form  solution  of  Eq.  (2)  can  be  derived  through  the 
QR-decomposition  of  A using  the  Gram-Schmidt  process.  In  the  following, 
we  denote  by  a*  and  qj  the  ith  column  vector  of  the  matrices  A and  Q, 
respectively,  and  by  the  ijth  element  of  the  upper  triangular  matrix  R. 
There  is  no  summation  over  repeated  indices.  The  general  solution  is 

x = W*b,  (3) 


where  W is  a do  x no  matrix  with  column  vectors  Wj  given  by 


with 


no 

Wj  = Cjjqj  'y  ^ Cj/:  ckk  qk, 
k=i+ 1 

(j<i 
3-i  + y ' Cji  B.j  J . 

3= 1 

The  geometrical  quantities  Cij  are  defined  as 
„-i 

j-i 

^ij  = ( Cii  I'ij  "b  y ' Cjfc  Ckf • Tkj  J for  J 1 ^ TIq, 

k=i+ 1 


where 


rij 


— ( aj  • &j  - yr  rki  rkj  J for  j >i<  n0. 
rii  \ k=\  ) 


(4) 

(5) 

(6) 

(7) 

(8) 
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The  general  closed-form  solution  allows  the  reconstruction  of  derivatives 
to  an  arbitrarily  high  order  of  accuracy  on  unstructured  grids. 

3.  Least-Squares  Filtering 

The  least-squares  gradient-reconstruction  method  can  be  turned  into  a fil- 
tering method  by  modifying  Eq.  (1),  so  that  4>o  is  no  longer  a point  value, 
but  represents  a filtered  value  (f>0, 


(t>i  = (f> o + (V0)o  • Aroi.  (9) 

The  effect  of  this  modification  is  that  the  filtered  value  (f)0  is  appended  to 
the  vector  of  unknowns  x.  The  resulting  system  of  equations  can  be  solved 
using  the  method  described  in  Section  2.  This  leads  to  an  expression  for 
the  filtered  value  in  the  form  of  a weighted  sum, 

do 

4>0  = (10) 

i— 1 

The  accuracy  of  the  filtering  operation  is  determined  by  the  order  of 
the  derivatives  included  in  the  gradient  reconstruction.  For  example,  linear 
gradient  reconstruction  leads  to  a second-order  accurate  expression  for  the 
filtered  value,  with  weights  given  by 

1 ( , r23  A . n27"23  — H3U22  A \ 

Woi  = -o-  1 1 Ayoi  H Axoi  • (11) 

r33  \ r< 22  rnr22  ) 

It  is  easily  verified  that  Eq.  (11)  leads  to  two  vanishing  moments, 


do 

^w0jAroi  = 0, 

2=1 

which  is  merely  a consequence  of  its  second-order  accuracy. 

The  spectral  behaviour  of  Eq.  (10)  can  be  improved  if  the  unfiltered 
value  4>o  is  also  included  in  the  stencil, 


do 

<£o  = Wo  + (l-a,0o)£"o^-  (12) 

2=1 

This  modification  does  not  degrade  the  accuracy  of  Eq.  (10). 

The  remainder  of  this  work  is  based  upon  Eq.  (12)  and  investigates  lin- 
ear and  quadratic  filter  functions.  For  both  functions,  the  stencil  is  extended 
to  include  an  additional  layer  of  vertices  beyond  the  nearest  neighbours.  In 
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the  interior,  this  gives  a stencil  of  18  vertices  for  woo  = 0 or  19  vertices  for 
woo  7"  0.  The  motivation  for  extending  the  support  is  twofold.  First,  hav- 
ing two  layers  of  vertices  allows  the  introduction  of  two  parameters  (3  and 
7,  which  can  be  used  to  weight  the  contributions  to  Eq.  (10)  of  the  inner 
and  outer  layers,  respectively.  Together  with  woo,  these  additional  degrees 
of  freedom  may  be  used  to  optimize  damping  of  high-wavenumber  compo- 
nents or  to  achieve  a specified  filter  width.  Second,  the  nearest-neighbour 
stencil  leads  to  a singular  matrix  A for  the  quadratic  filter  function  on 
uniform  grids,  thus  necessitating  the  use  of  additional  points  in  the  stencil. 

One  advantage  of  the  new  filtering  method  is  that  it  allows  reusing 
data  structures  and  geometric  weights  already  employed  to  compute  gra- 
dients in  the  flow-solution  method.  Furthermore,  it  is  easily  extended  to 
three  dimensions  and  does  not  require  special  treatment  at  boundaries  be- 
yond ensuring — as  for  interior  vertices — that  the  degree  of  a given  vertex 
is  greater  than  the  number  of  derivatives  reconstructed  at  that  vertex. 


4.  Determination  of  Filter  Width 

In  the  present  work,  the  filter  width  is  determined  from  the  polar  moment 
of  the  filter  transfer  function, 


it/ A tt/A 

JXy  = J J (kl  + ky')G(kx,ky)dkxdky,  (13) 

o o 


where  k = { kx , kyY  is  the  wave-number  vector,  G(kx,  ky)  is  the  filter  trans- 
fer function,  and  A is  the  grid  spacing.  In  one  dimension,  this  reduces  to 
the  second  moment  of  the  filter  transfer  function,  whose  use  in  determining 
the  filter  width  was  originally  suggested  by  Lund  (1997).  In  this  section, 
we  assume  the  grid  spacing  to  be  uniform. 

The  ratio  of  the  filter  width  A/  to  the  grid  spacing  may  be  computed 
from  the  relation 


8 JXyA4  J 


(14) 


The  constants  in  Eq.  (14)  were  chosen  such  that  it  gives  the  correct  width 
for  the  Fourier  cut-off  filter. 

Figure  1(a)  depicts  the  transfer  function  for  the  linear  filter  function 
with  woo  — 1/5  and  /3  = 7 = 1.  While  high  wave-numbers  are  damped 
well,  the  transfer  function  deviates  quickly  from  unity  for  low  to  moderate 
wave-numbers,  and,  as  such,  is  a poor  representation  of  the  Fourier  cut-off 
filter.  Numerical  evaluation  of  Eq.  (13)  gives  a = 1.40. 
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The  transfer  function  for  the  quadratic  filter  function  with  woo  = 1/2 
and  (3  = 7 = 1 is  shown  in  Fig.  1(b).  Compared  to  the  linear  filter  func- 
tion, the  quadratic  filter  function  is  a good  approximation  to  the  Fourier 
cut-off  filter  for  low  to  moderate  wave-numbers.  For  higher  wave-numbers, 
preferred  directions  can  be  discerned  which  are  aligned  with  the  edges  in 
the  grid.  For  the  quadratic  filter  function,  Eq.  (13)  gives  a = 1.12. 


Figure  1.  Transfer  functions  for  filter  functions  on  uniform  grids  and  /?  = 7 = 1.  (a) 
Linear  filter  function  with  u>oo  = 1/5  and  (b)  quadratic  filter  function  with  cuoo  = 1/2. 
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5.  Commutation  Error 

Marsden  et  al.  (2000)  proved  that  a filter  with  p — 1 vanishing  moments 
is  needed  to  achieve  a commutation  error  of  order  p for  smoothly  varying 
filter  widths  in  one  dimension.  An  equivalent  statement  is  that  the  filter 
must  be  accurate  to  order  p to  achieve  a commutation  error  of  order  p 
for  smoothly  varying  filter  widths.  For  grids  with  arbitrarily  varying  filter 
widths,  the  commutation  error  will  drop  below  p.  It  is  thus  necessary  to 
construct  filter  operators  of  order  p + 1 to  obtain  commutation  errors  of  or- 
der p on  arbitrary  unstructured  grids.  We  are  interested  in  this  general  case 
since  unstructured  grids  rarely  satisfy  smoothness  constraints.  To  obtain  a 
second-order  commutation  error,  we  thus  require  quadratic  filtering. 

The  order  of  the  commutation  error  is  computed  by  carrying  out  a 
grid-refinement  study  using  an  analytic  function  for  the  unfiltered  field. 
Five  uniform  triangular  grids  were  generated  for  a hexagonal  domain,  con- 
taining 271,  1141,  4681,  18961,  and  76321  vertices.  The  interior  vertices 
were  subsequently  distorted  by  random  amounts  of  a given  fraction  of  the 
grid  spacing  in  both  coordinate  directions.  In  the  results  presented  below, 
this  fraction  was  taken  to  be  0.35.  The  distorted  grid  with  1141  vertices  is 
shown  in  Fig.  2. 


Figure  2.  Distorted  grid  with  1141  vertices.  Inset  shows  detail  of  distorted  grid. 
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The  commutation  error  is  defined  in  terms  of  the  discrete  divergence, 


/ 6u  6v  \ / 5u  foA 

y 8x  8y  J \8x  + SyJ  ’ 


(15) 


where  8(-)/8x  represents  the  discrete  gradient  operator  in  the  ^-coordinate 
direction.  The  function  chosen  in  the  present  study  is 


J u 1 _ J cos(7nc)  sm(ny)  1 
\ v j \ sin(7nr)  cos(7ry)  J ' 


(16) 


In  computing  the  commutation-error  norms,  only  the  vertices  were  in- 
cluded whose  gradients  or  filtered  values  were  not  affected  by  boundary 
effects.  Effects  of  one-sided  stencils  arising  from  the  presence  of  boundaries 
will  be  studied  in  future  work.  It  was  verified  that  commutation  errors  were 
identically  zero  for  arbitrary  functions  on  uniform  grids. 

The  variation  of  the  /^-norrns  of  the  commutation  error  with  grid  re- 
finement is  shown  in  Fig.  3.  Note  that  the  commutation  error  is  about 
an  order  of  magnitude  smaller  than  the  truncation  error  of  the  divergence 
operator.  The  order  of  accuracy  of  the  filtering  operator,  the  divergence 
operator,  and  the  order  of  the  commutation  error  were  computed  from  a 
linear  least-squares  curve  fit  for  the  finest  four  grids.  The  slopes  were  de- 
termined to  be  3.05,  1.96,  and  2.17,  respectively.  It  is  thus  verified  that  the 
commutation  error  obtained  with  the  new  quadratic  filtering  method  on  a 
randomly  distorted  unstructured  grid  is  of  second  order. 

The  ultimate  test  for  commutation  will  be  to  specify  a uniform  filter 
width  on  a randomly  distorted  grid  using  the  parameters  u>oo,  /?,  and  7 and 
to  check  for  zero  commutation  error.  This  is  an  objective  of  future  work. 


6.  Conclusions 

A new  filtering  method  for  unstructured  grids  was  presented.  Closed-form 
expressions  were  given  which  allow  the  construction  of  filtering  operators  of 
arbitrarily  high  order.  The  new  filtering  method  is  easily  constructed,  does 
not  require  special  treatment  at  boundaries,  and  allows  reusing  data  struc- 
tures and  geometric  terms  needed  by  the  flow-solution  method  without  fil- 
tering. Linear  and  quadratic  filter  functions  were  studied.  A grid-refinement 
study  on  randomly  distorted  grids  demonstrated  a commutation  error  of 
second  order. 
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Figure  3.  Variation  of  Z/2-norm  of  errors  with  grid  refinement.  Solid  lines  represent 
linear  curve  fits  to  data,  and  N denotes  the  number  of  vertices. 
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